Thermodynamic properties of the Yb 2 Ti 2 07 pyrochlore as a function of temperature 
and magnetic field: validation of a quantum spin ice exchange Hamiltonian 
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The thermodynamic properties of the pyrochlore Yt)2Ti2 07 material are calculated using the 
numerical linked-cluster (NLC) calculation method for an effective anisotropic-exchange spin-1/2 
Hamiltonian with parameters recently determined by fitting the neutron scattering spin wave data 
obtained at high magnetic field h. Magnetization, M(T, h), as a function of temperature T and for 
different magnetic fields h applied along the three high symmetry directions [100], [110] and [111], 
are compared with experimental measurements on the material for temperature T > 1.8 K. The 
excellent agreement between experimentally measured and calculated M(T, h) over the entire tem- 
perature and magnetic field range considered provides strong quantitative validation of the effective 
Hamiltonian. It also confirms that fitting the high-field neutron spin wave spectra in the polarized 
paramagnetic state is an excellent method for determining the microscopic exchange constants of 
rare-earth insulating magnets that are described by an effective spin-1/2 Hamiltonian. Finally, we 
present results which demonstrate that a recent analysis of the polarized neutron scattering inten- 
sity of Yb2Ti2C>7 using a random phase approximation (RPA) method [Chang et al, Nat. Comm. 
3, 992 (2012)] does not provide a good description of M(T, h) for T < 10 K, that is in the entire 
temperature regime where correlations become non- negligible. 

PACS numbers: 75.10.Jm, 75.40.Cx,75.47.Lx,75.30.Et 



I. INTRODUCTION 



Magnetic rare-earth pyrochlore oxides and effective 
spin-1/2 quantum dynamics 



Quantum spin liquids are magnetic systems in which 
large quantum mechanical zero-point spin fluctuations 
prevent the development of long-range order down to 
absolute zero temperature. The search for real mate- 
rials that display this phenomenology in two and three 
dimensions is a very active research topic in the field 
of condensed matter physics X The interest in spin liq- 
uids stems from the expectation that they may host non- 
trivial quantum entanglement, topological order, as well 
as emergent fractionalizcd and dcconfined low-energy ex- 
citations. Spin liquids have also been conjectured as 
progenitors of unconventional superconductivity at "high 
temperature" . Highly frustrated magnets are particu- 
larly apt at exhibiting large quantum spin fluctuations. 
These magnets are realized in systems which consist of 
localized magnetic moments (spins) that reside on two 
and three dimensional lattices of corner-sharing triangles 
or tetrahedra and which interact with effective antifcrro- 
magnetic nearest-neighbor coupling.— Such highly frus- 
trated magnets display an exponentially large number 
of classical ground states^ This allows for quantum me- 



chanical effects to be tremendously magnified compared 
to magnetic systems with conventional long-range mag- 
netic order ji 

The pyrochlore oxides, of chemical formula 
RE2M207^ count a multitude of magnetic members. 
In RE2M2O7, RE is a magnetic trivalent (Lanthanidc, 
4f) rare-earth ion (Gd 3+ , Tb 3+ , Dy 3+ , Ho 3 +, Er 3 +, 
Yb 3+ ) or a non-magnetic ion (Y 3+ , Lu 3+ ), while M 
is tetravalent, typically a transition metal ion, which 
can be either mag netic (Mo 4+ , Mn 4+ ) or not (Ti 4 +, 
Sn 4+ , Zr 4+ , Gc 4+ ). Of relevance to the above discussion 
is the fact that the RE and M ions reside on two 
distinct and interpenetrating three-dimensional lattices 
of corner-sharing tetrahedra. As a result, high geometric 
magnetic frustration ensues whenever the RE-RE or M- 
M interaction is effectively antiferromagnetic^. Among 
the RE2M2O7 family, the magnetic rare-earth oxides, in 
which RE is magnetic and M is not, have been rather 
extensively studied^ These have revealed a number of 
fascinating phenomena such as long-range order induced 
by order-by-disorder^ 7 - multiple-fc long-range ordered 
phased spin liquid behavior— and spin ice physics^ - — 
the latter having attracted much interest i^i 5 . : 13 ) 14 The 
spin ice state displays a residual low-temperature 
magnetic entropy close to that found for common water 
ice ; 12 : 15 : 16 hence the name spin ice. Most recently, a 
renewed flurry of theoretical and experimental efforts 
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have been directed at the study of spin ices for they have 
been argued to display an emergent Coulomb phased 
accompanied at low energies by deconfincd fractionalized 
magnetic charge excitations, or "monopoles" 1 14 i 17 i 18 

Notwithstanding the rich physics that RE2M2O7 ma- 
terials display^ they have, until very recently^ - — not 
attracted that much interest from the community of the- 
orists and experimentalists searching for quantum spin 
liquidsii Quantum fluctuations are expected to be, rel- 
atively speaking, more significant the smaller the spin 
quantum number S. On the other hand, because of 
the intrinsic large spin-orbit coupling at play in the Lan- 
thanide 4f series, the total angular momentum J = L+S, 
and not S, is a good quantum number with J being 
typically large across the whole Lanthanide series (e.g. 
J = 8 and J = S = 7/2 for free Ho 3 + and Gd 3 + ions, 
respectively). This situation would seem unpromising to 
those seeking magnetic materials with potentially large 
quantum mechanical spin fluctuations and exotic quan- 
tum states of matter. However, such a perspective is 
perhaps too pessimistic as we now discuss. 

In insulating magnetic compounds with 4f elements, 
the various inter-ion interactions, i/mt (exchange, su- 
perexchange, virtual phonon exchange and magneto- 
static dipolc-dipolc coupling) are typically weak com- 
pared to the single-ion crystal-field interactions defined 
by a crystal- field Hamiltonian H c f. As a first approxi- 
mation, one thus often proceeds by determining the en- 
ergy level spectrum of H c f for a fixed J manifold^ The 
point group symmetry of the crystal dictates the allowed 
symmetry properties of H c {. These symmetries deter- 
mine, in return, the spectral decomposition of the crystal- 
field states that derive from the original 2 J + I degener- 
ate 2S+1 Lj electronic ground state of the otherwise free 
RE 3+ ion. In the simplest case, the resulting crystal 
field ground state is a magnetic doublet with wavefunc- 
tions \4> + ) and \ip~) that have = J2 m , Cm,\J, m j) 
for spectral decomposition, \tp^) contains all the \ J,mj) 
spectral components that transform similarly according 
to the point group symmetry operations. Consequently, 
there may or not be nonzero (-0 ± |iJi nt |'0 :F ) matrix ele- 
ments. This depends on the specific nature of the inter- 
ion couplings -ffint(Jj U ), a function of the components 
Jf (u = x, y, z) of angular momentum Ji of ion 
as well as the specific ion-dependent spectral decomposi- 
tion of |V' ± )- It is the nonzero (V' ± |-ffint|V' =F ) matrix el- 
ements^ which determine whether significant quantum 
dynamics exist within the low-energy sector. Most im- 
portantly, quantum dynamics need not be ruled out de- 
spite the large J of the isolated RE 3+ since the crystal 
field Hamiltonian H c { entangles a superposition of the 
mj) eigenstates of J z 2£ As a result, H int , by virtue of 
its lack of commutation with H c f , can, in principle, have 
nonzero matrix elements between \ip + ) and \ip~) and in- 
duce quantum dynamics within the low-energy Hilbert 
space spanned by f]f VV^I^D- 

In 4f ions with an even number of electrons (i.e. non- 
Kramers ion) such as Pr, Tb and Ho, time-reversal 



symmetry imposes that (?/» ± | J ± |-0 :F ) = 0<2i Conse- 
quently, in presence of solely bilinear interactions of the 
form Kl l J v (r i j)Jl l Jj with anisotropic Kf? couplings, non- 
Kramers ions would display no quantum dynamics at 
low-energy and behave as effective classical Ising spins 
S=I/2, as in the well-studied LiHoF4 dipolar Ising sys- 
tem^ In the (Pr,Tb,Ho) 2 (Ti,Sn,Zr) 2 7 materials^ in- 
teractions beyond bilinear ones or consideration of the 
excited crystal field states are necessary to cause quan- 
tum dynamics in the low-energy sector j^2t— In that con- 
text, multipolar interactions in Pr2(Sn,Zr,Ir)207 com- 
pounds^ and virtual crystal field excitations in the 
Tb 2 (Ti,Sn) 2 07 materials have been discussed i 27 i 28 Con- 
versely, odd-electron (Kramers) ion systems (e.g. Gd, 
Dy, Er, Yb) are in principle symmetry-allowed to have 
nonzero (■0 ± | J ± |-0 =F ) matrix elements. The famous 
Dy 2 Ti 2 07 spin ice compound, in which the Dy 3+ ions 
are Kramers ions and for which the excited crystal 
field states lie at ~ 300 K above the ground doublet, 
has been shown to be well-described by a dipolar spin 
ice mode l 30 ! 31 with classical Ising spins4 This success 
very likely signals rather negligible interactions among 
the Jf components beyond bilinear ones, concomitantly 
with the specific spectral decomposition of |^/» ± ) for 
Dy 3+ in Dy 2 Ti 2 7 4 2 .^ On the other hand, Er 2 Ti 2 7 
and Yb 2 Ti 2 07 have been known for some time^ ' 33 ' 34 to 
have predominant "transverse" ("0^1 J ± |'0 =F ) matrix ele- 
ments^ along with a non-negligible (ip ± \J z \ip ± ) "lon- 
gitudinal" {g zz tensor) component. Yet, the two com- 
pounds display quite different behaviors. Er 2 Ti 2 07 has 
overall antifcrromagnctic interactions and develops long- 
range order at 1.2 K with zero propagation vector q or d 
and zero magnetic moment per tetrahedron^ that is in- 
duced by order- by-disorder iL In contrast, Yb 2 Ti 2 07 has 
overall ferromagnetic interactions and exhibits a phase 
transition at T c ~ 0.24 However, the nature of the 

long-range order below T c remains disputed 3 ^— and the 
high-sensitivity of the properties in the low-temperature 
regime (T < 300 mK) on sample quality is just beginning 
to be understood.— ~— It is here, within the Yb 2 M 2 07 
family with overall ferromagnetic interactions and signifi- 
cant (i/' ± | J^lip^) matrix elements ; 20 ' 24 that the potential 
for an exotic class of quantum spin liquid arises - a pos- 
sibility that may have been casually dismissed from the 
naive perspective of "there should be negligible quantum 
effects in such large J systems" . 



Quantum spin ice and Yb2Ti2C>7 

The Dy 2 (Ti,Sn,Ge) 2 7 and Ho 2 (Ti,Sn,Gc) 2 7 mate- 
rials are classical Ising systems which may be viewed in 
their spin ice regime as collective paramagnets* 3 - or, em- 
ploying a more contemporary terminology, classical spin 
liquids.— As discussed above, highly frustrated magnetic 
systems, by virtue of their low-propensity to develop clas- 
sical long-range order, are attractive candidates to search 
for quantum spin liquid behavior. Spin ice, reinterpreted 
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as a classical spin liquid) 1 - may thus be viewed as a natural 
setting to explore how the addition of quantum dynam- 
ics may give rise to a quantum spin liquid state. Such 
a topic was originally explored by Hermele, Fisher and 
Balents^ and Castro-Neto, Pujol and Fradkini 5 - a few 
years ago in the context of minimal theoretical toy mod- 
els. These two groups argued in their respective paper 
that the addition of quantum dynamics within a par- 
ent (constrained) classical spin ice manifold may pro- 
mote the system to a U(l) quantum spin liquid. Such 
a state would be dcscribablc by a quantum field theory 
analogous to that of quantum electrodynamics (QED) in 
3+1 dimensions. As a consequence, this U(l) spin liq- 
uid state would display "electric" and "magnetic" charge 
excitations and an accompanying gauge boson, or "ar- 
tificial photon" i 44 i 45 Numerical simulations have, over 
the past few years, suggested that such phenomenology 
may indeed be at play in various lattice toy models . 46 ' 47 
Such a U(l) quantum spin liquid state, which may be 
referred to as "quantum spin ice" j2! has been suggested 
to explain some of the properties of real materials such 
as TbaTiaOT^s an d Pr 2 (Sn,Zr) 2 7 .-S2, Most recently, it 
has been proposed that Yb-based Yb2M2C>7 pyrochlorc 
oxide s 20 ' 21 may offer an exquisite class of systems to in- 
vestigate the possible existence of a quantum spin ice 
state and where the complexities of virtual crystal field 
fluctuation s 27 ' 28 and magnetoelastic couplin g 48 ' 49 that 
complicate the Tb2(Ti,Sn)2C>7 compounds are avoided. 
In a theory paper building on the original work of Her- 
mele et al.^ Savary and Balents have recently put for- 
ward a mean-field lattice gauge theory which identifies 
a number of possible phases, the most exotic ones being 
the aforementioned U(l) spin liquid as well as a novel 
Coulomb ferromagnetic state.— This approach has been 
further extended to non-Kramers ion systems^ Finally, 
the question of how to probe the emergent photon in 
quantum spin ice via inelastic neutron scattering mea- 
surements has very recently been discussed.— 

One particularly interesting aspect of the search for 
quantum spin liquids in the Yb-based pyrochlores, and 
perhaps Pr-based pyrochlores as well , 22 ' 29 ' 51 is that the 
microscopic Hamiltonian is parametrized by a handful of 
independent anisotropic exchange parameters { J e } (four 
for Yb 3+ and three for Pr 3+ ) 2H£2rJ& between effective 
spin-1/2 degrees of freedom on each pyrochlore lattice 
site. Furthermore, it may be that the long-range dipolar 
interactions, so important in classical Dy2(Ti,Sn)207 and 
Ho2(Ti,Sn) 2 07 spin ices with the large magnetic Dy 3+ 
and Ho 3+ magnetic moment,- ' 11 ' 30 ' 31 ' 55 ' 56 can perhaps 
be neglected as a first approximation. To make definite 
progress at this time, a determination of the effective 
anisotropic exchange from experiments and theory using 
controlled approximations is required in order to posi- 
tion a candidate spin ice material in the phase diagram 
of Ref. [2l[, or the corresponding phase diagram rele- 
vant to non-Kramers ions^ In the case of Yb 2 Ti 2 07, 
the determination of those interactions from a series of 
experiment s 20 ' 39 ' 52 ' 57 ' 58 has led to different {J e } sets with 



very different numerical values. Perhaps most notewor- 
thy, a determination of those parameters from a fit to spin 
waves in strong magnetic field measured via inelastic neu- 
tron scattering 2 ^ give values significantly different than 
those obtained by fitting the zero-field diffuse paramag- 
netic neutron scattering using a random phase approxi- 
mation (RPA) method 5 ^ at a temperature T ~ 1.4 K or a 
subsequent polarized neutron scattering version, also us- 
ing RPA as the fitting procedure, but now very near T C M 
Encouragingly, however, recent numerical linked-cluster 
(NLC) calculations^ 9 - have found the zero-field magnetic 
specific heat data of Yb2Ti2C>7 to be well described above 
0.7 K by the { J e } set obtained from the inelastic neutron 
scattering^ - but not by the other sets determined from 
RPA fits to diffuse neutron scattering ! 39 ' 52 

It is clearly desirable to fit bulk measurements to deter- 
mine the {J e } parameters in cases where inelastic neu- 
tron scattering studies are impractical, and to corrob- 
orate such neutron studies and understand thermody- 
namic and bulk magnetic properties in a common frame- 
work with neutron studies whenever possible. Yet, the 
scarcity of controlled numerical methods readily available 
to calculate the thermodynamic properties of frustrated 
three-dimensional anisotropic quantum spin systems in 
a temperature regime where non-trivial correlations de- 
velop do not make this task straightforward. With the 
seeming success of previous NLC calculations applied to 
Yb2Ti2 07)£ 9 - we are thus naturally led to ask: 

Can one convincingly demonstrate that NLC does pro- 
vide a controlled method to describe bulk data for a three- 
dimensional frustrated quantum (spin ice) system as it 
progressively enters its low-temperature strongly corre- 
lated regime? 

In this paper, we address this question by extending 
the work of Ref. [59| by computing the thermodynamic 
properties of Yb2Ti207 in nonzero magnetic field using 
NLC and by comparing the results of such calculations 
with measurements above 2 K on single crystals. We 
show, through a comparison with NLC, that the effec- 
tive Hamiltonian with its anisotropic exchange parame- 
ters {J e } previously determined by fitting spin waves in 
the field polarized paramagnetic stated - describes, with 
no adjustment, the temperature T and magnetic field h 
dependence of the magnetization, M(T, h), of Yb 2 Ti207. 
As a consequence, we demonstrate simultaneously the 
usefulness of the NLC method and further validate^ 9 - the 
quantitative merit of the effective spin Hamiltonian of 
Ref. [U to describe Yb 2 Ti 2 7 . 

The rest of the paper is organized as follows. In the 
next section, we describe the effective spin-1/2 model for 
Yb 2 Ti 2 7 along with the NLC method. In Section llHl 
we discuss the details of the experimental method em- 
ployed. The NLC results are presented in Section ITVl 
while Section [V] provides a comparison between experi- 
ment and theory. A brief discussion in Section I VII con- 
cludes the paper. 
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II. MODEL & COMPUTATIONAL METHOD 

Symmetry considerations imply that nearest-neighbor 
bilinear exchange constants on the pyrochlore lattice can 
be parametrized in terms of 4 distinct exchange param- 
eters and two g-factors. In the local basis, this model is 
defined by the Hamiltonia n 20 ' 21 ' 51 " ^ 

%qsi = y^,{JzzS?Sj ~ J±(S+S~ + S~S+) 

+J±±htjS+S+ + jfjS^Sr] 
+J,±[(Sf(C«S+ + (i 3 S~) + i o j}}. (1) 

Several notation conventions are possible for 
Hqsi2&22t— Here, we adopt the one used in Ref. 
In Eq. ([1]), refers to nearest-neighbor sites of the 

pyrochlore lattice, jij is a 4 x 4 complex unimodular 
matrix, and C = — 7*i 20 ' 21 The z quantization axis is 
along the local [111] direction^ and ± refers to the 
two orthogonal local directions. The g tensor takes 
values g zz along the local [111] cubic direction and g xy 
perpendicular to it. In the presence of an applied exter- 
nal magnetic field h 7 an additional Zeeman interaction 
Hz = — h ■ y ■ SfiB is added to 'Hqsi, giving a total spin 
Hamiltonian H. = Hqsi + Tiz- Here is the g-tensor 
with eigenvalues (g xy ,g xyi g zz )M^- 

As already stated above, there have been several 
experimental attempts to determine these {J c } effec- 
tive exchange parameters leading to widely different re- 
sults i^^^iS 2 -! 5 -^ Among them, Ross et alM used in- 
elastic neutron scattering (INS) data obtained from mea- 
surements in high magnetic fields (i.e. in the polar- 
ized paramagnetic state) to determine the { J e } couplings 
for Yb 2 Ti 2 7 : J zz = 0.17 ± 0.04, J± = 0.05 ± 0.01, 
J±± = 0.05 ± 0.01, and J z ± = -0.14 ± 0.01, all in meV 
and with g zz = 1.80 and g xy = 4.32. In a recent study^ 
it was found that the zero-field specific heat and entropy 
deduced from the data of Blote et al^ are well described 
only by the exchange parameters obtained by Ross et 
al2& A key goal of the present paper is to investigate fur- 
ther the validity of the Hamiltonian (flj as a description 
of the bulk properties of YbgTiaO? with the {J e } ex- 
change couplings from Ref. |2(j] ■ We do so by performing 
calculations of the thermodynamic properties of model 
(1) as a function of temperature T and magnetic field h 
and compare the numerical results with those obtained 
from experimental measurements. We find an excellent 
agreement without any adjustment of the parameters de- 
termined by Ross et aZ.jSP. hence providing a strong vali- 
dation to the {Je] parameters determined by INS. As a 
corollary, our work confirms that fitting the INS spectra 
at high magnetic field is an excellent method to deter- 
mine the exchange constants for pyrochlore oxides with 
well isolated magnetic ground doublets and which are 
well described by an effective spin-half modcl^ Further 
implications of this agreement will be discussed in the 
concluding Section IVT1 



To calculate the thermodynamic properties of model 
([T]), we turn to Numerical Linkcd-Cluster (NLC) expan- 
sions^ In this method, an extensive property P (such 
as heat capacity or magnetization) of a thermodynamic 
system is calculated as a sum over contributions from 
different clusters embedded in the lattice: 

P/N = ^L(c) xW(c). (2) 

c 

Here L(c) is the count of the cluster per lattice site, de- 
fined as the number of ways to embed the cluster. W(c) 
is the weight of the cluster which is obtained by calculat- 
ing the property for a given cluster and subtracting the 
weight of all its subclusters. 

W(c)=P(c)-J2w(s). (3) 

Here, the sum runs over all subclusters of the cluster c. 
A subcluster is defined as any cluster smaller than c that 
can be embedded in cluster c. This scheme can be used 
to develop power series expansions, such as high temper- 
ature series expansions in powers of inverse temperature 
/3, or, some coupling constant expansion (such as expan- 
sion in inverse field-strength). It can also be used to 
numerically calculate properties for a given value of tem- 
perature and coupling constants, as we do in the present 
work. 

The pyrochlore lattice is a lattice of corner-sharing 
tetrahedra and it proves useful to develop NLC in terms 
of clusters consisting of complete tetrahedra. We have 
calculated temperature and field-dependent properties 
up to 4 th order, that is, including contributions from all 
clusters made of up to 4 tetrahedra. These NLC cal- 
culations are numerically exact in two limits. They are 
so at high temperatures, since corrections to 4 th order 
NLC is of order /3 6 in the high temperature series ex- 
pansion for hiZ. The NLC calculations are also exact 
at high magnetic field h at all temperatures since cor- 
rections to 4 th order NLC for hiZ is of order (J/h) 5 at 
T = with exponentially small corrections exp(— ch/T) 
at finite temperatures (c is some h and T independent 
constant). The parameter region where NLC begins to 
lose accuracy is when the temperature and the applied 
field (Zeeman energy) are both smaller than the exchange 
constants. Thus, as long as either the high temperature 
expansion or the high-field expansion converges, the NLC 
calculations should be highly accurate. 

The reason for developing NLC in terms of complete 
tetrahedra comes from the fact that, for spin ice systems, 
the 'ice rule' constraints at the origin of spin ice physics 
are local to tetrahedra. Thus, having clusters with only 
parts of a tetrahedron would cause wild oscillations in 
the calculated properties as the constraints cannot be 
satisfied for such a cluster. Having only clusters with 
full tetrahedra allows the system to always satisfy the 
ice-rule constraint. For example, in the case of the clas- 
sical nearest-neighbor spin-ice exchange model, such a 
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tctrahedra-based NLC was found to be highly accurate^ 
In fact, the first order NLC, which uses a single tetra- 
hedron and is equivalent to Pauling's approximation for 
the entropy of ice, gives quantities at all temperatures (in 
zero external field) that are accurate to a few percent^ 
However, in the quantum spin ice problem, NLC must 
break down at low temperatures due to the development 
of either long-range correlations or long-range entangle- 
ment. We have found (see results below and in Ref. [591 ]) 
that the first signature of such a breakdown in the low 
temperature and low field region is an alternation of the 
thermodynamic property P considered. This reflects the 
fact that, as system sizes increase, the thermodynamic 
properties must approach their infinite size values in some 
way. Now, imagine that, at some order, the numerically 
calculated property is a bit too large, giving rise to a large 
cluster weight W. Subcluster subtraction then causes the 
weight in the subsequent order to be too small. This, in 
return, causes an alternation with the NLC order consid- 
ered in the property obtained by restricting the sum to 
some order. Such an alternation is handled well by using 
the Eulcr summation method^ which ensures that an al- 
ternating piece is completely eliminated from the partial 
sums at all orders. Thus, the Euler resummed properties 
are only missing the longer range correlations which are 
necessarily absent in finite order NLC. For Yt>2Ti207, we 
have found that in zero magnetic field^ the thermody- 
namic properties converge down to 2 K without the Eulcr 
summation and down to about I K with Euler summa- 
tion (the largest exchange constant for Yb2Ti20y, J zz in 
Eq. Q]), is approximately 2 K— ). 



Computational details 

NLC using the tctrahcdral basis requires exact diago- 
nalization of increasingly larger tetrahedral clusters. Us- 
ing a single Intel® Core i7-920 processor and freely- 
available LAPACK linear algebra routines, diagonaliza- 
tions for clusters of one tetrahedron (four sites) and two 
tctrahcdra (seven sites) could be done in less than a sec- 
ond, while the three-tetrahedron (10-site) cluster still re- 
quired less than 10 seconds. Computing only the spec- 
trum for a single four-tetrahedron (13-site) cluster re- 
quired approximately 25 minutes of cpu time and 2.1 
GB of memory, or slightly over twice the memory re- 
quired to store the full Hamiltonian matrix (2 26 complex 
numbers). Generating the full set of eigenstates required 
between 4 and 8 GB of memory. With the computer re- 
sources available to us, full exact diagonalization of larger 
systems, and thus higher orders of NLC expansion, were 
prohibited by memory requirements. 

A list of different topological clusters with complete 
tetrahedra is shown in Fig. [TJ A single site must be 
treated as a cluster as well, but is not shown. Because 
the local quantization z axis varies from sublattice to 
sublattice in the RE2M2O7 pyrochlorc oxidcs^2^ in the 
presence of a magnetic field, each cluster with different 
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FIG. 1: (color online). List of topological clusters with com- 
plete tetrahedra. In an applied field, each cluster must be 
treated as several distinct clusters as explained in the text. 
The integer n refers to the order at which the cluster arises 
and L gives the cluster count in zero-field. 



types of sublattice site must be treated separately. For 
example, each of the four single-site graphs (in NLC order 
n = 0) must be treated separately for a general field di- 
rection. A single tetrahedron graph remains unique even 
in a field. But, the graph with two tetrahedra joined at 
a point (NLC order n = 2), is really 4 distinct graphs 
that depends on the type of site (i.e. sublattice) that is 
shared between the two tetrahedra. Similarly, all higher 
order graphs must be split into several distinct graphs to 
complete the calculation. 

Calculation of observables 

We use the eigenvalue spectrum {E a } of% = "Hqsi + 
Hz to compute the partition function Z: 

Z = tr(e-W) = J2^ Ea - ( 4 ) 

a 

The expectation value for the internal energy is computed 
from the formula 

{E)=tx{He- m )/Z. (5) 
This quantity is used, in turn, to compute the entropy as 
(S)/k B =/?(£)+ log (6) 
to which the heat capacity is related as 

(C)=4(S). (7) 
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The average magnetization per site M(T, h) is calcu- 
lated as a derivative of the free energy with respect to 
the applied magnetic field h: 

M{T,h)=p- 1 ^\nZ. (8) 

In practice, an approximation of this derivative at a given 
field value requires two separate evaluations of the free 
energy separated by a small change of field strength, in 
this case 1CT 6 T. 

We believe in the correctness of the results from our 
calculations of the field-dependent properties for the fol- 
lowing three reasons. First, the same results were ob- 
tained from two independently written computer pro- 
grams. Second, we checked that in the limit h — > 0, the 
zero-field results 5 - 9 - were reproduced. Third, we checked 
that at high temperatures the weights of the clusters were 
found to scale with the expected powers of inverse tem- 
perature /3. These powers can be deduced from consid- 
erations of a high temperature expansion of the quan- 
tity of interest. Consider, for example, an expansion for 
InZ. In such an expansion either at least two powers of 
(3=1/ (/sbP) arise for each tetrahedron or at least one 
power of (3 arises for each tetrahedron together with two 
powers of j3 from placing the Zeeman T~Lz field-term on 
sites at the outside perimeter of the cluster, as needed to 
give a non-zero contribution to the trace in Eq. 2) Thus, 
the weight of the 4-tetrahedra cluster is of order /3 8 in zero 
field and of order /3 6 in non-zero field. These latter checks 
are non-trivial and demonstrate that all subgraphs have 
been properly subtracted within the NLC procedure. 

Euler summation 

NLC generates a sequence of property estimates {P n } 
with increasing order n, where P n = X^j=i The con- 
vergence of such a sequence can be improved by Eu- 
ler summation! 59 ! 62 In general, given alternating terms 
Si = (—iyui, the precise infinite-size lattice property 
Poo (£) is approached by the sum (with n even) 

°° (-l) s 

Uq - Ul + U 2 - . . . - U n -1 + ^2 2S+1 [A S tt„], (9) 

where A is the forward difference operator 
A°u n = u n , 

Au n = U n +1 - U n , 
A 2 U n = U n+ 2 - 2u n+ i + U n , 

A z u n = u n+3 - 3u n+2 + 3u„ + i - u n , (10) 

Usually, a small number of terms are computed directly, 
and the Euler transformation, P nj E defined below, is ap- 
plied to the rest of the series. In our case, where direct 
terms were available up to fourth order, we began the 
transformation after the second order, so that the third 



and fourth order Euler-transformcd property estimates, 
P3E and P^e respectively, are given by 

P3,E = <5o + Si + S-2 + 2^*3' 

P4,E = P 3 ,E + ^^. (11) 

III. EXPERIMENTAL METHODS 

As we are aiming to establish a close connection be- 
tween experiment and theory for the temperature and 
field dependence of the magnetization M{T, h), it is nec- 
essary to perform adequate demagnetization corrections 
and this, in turn, necessitates measurements on single 
crystals. Magnetization experiments were thus carried 
out on three single crystal pieces. The samples were 
cut from two large single crystals grown by the opti- 
cal floating zone method, using a growth procedure sim- 
ilar to that previously employed for growing single crys- 
tals of Tb2Ti 2 One crystal, which provided samples 
aligned along [110] and [100], was grown at a rate of 
6mm/h in oxygen pressure of 4atm. The second crystal, 
from with a piece aligned along [111] was cut, was grown 
at 5mm/h in 2atm of oxygen. The samples were aligned 
using X-ray diffraction to within 2 degrees of of each 
of the three high-symmetry directions; [100], [110], and 
[111]. Two single crystal samples were cut into rectan- 
gular prisms ("needles") measuring 0.74 mm x 0.74 mm 
x 3.74 mm and 0.69 mm x 0.65 mm x 2.67 mm, respec- 
tively, with [111] and [100] directions oriented along their 
long axes. The applied magnetic field was also oriented 
parallel to the long axis. The third single crystal piece 
was cut and polished into a rounded triangular shape 
which we approximate, for the purposes of demagneti- 
zation corrections, as an ellipsoid with major axes a = 
1.85 mm, b = 1.5 mm and c = 0.8 mm. The a direction 
made an angle of 22° with respect to both [110] and the 
applied field direction within the ab plane. The magne- 
tization of the same [110] sample was also reported in 
Ross et al.r^ but without a demagnetization correction 
applied. The data was collected with a Quantum Design 
MPMS instrument, which uses a SQUID magnetometer 
to measure the DC magnetization in magnetic fields up to 
5 T and temperatures as low as 1.8 K. The magnetization 
data for the [100] and [111] samples were corrected for de- 
magnetization effects using an approximate formula for 
the demagnetization field in rectangular prisms^ 5 - The 
data for the [110] sample was corrected by approximat- 
ing it as a very flat ellipsoid, with the field direction 22° 
from the (long) a direction. In general, the susceptibility, 
M(T,h)/h, can be corrected to account for the demag- 
netization field using the following formula, in SI units, 

l/ X =l/XA-N (12) 

where x is the actual susceptibility that we aim to com- 
pare with the results of NLC calculations, xa is the ap- 
parent susceptibility (given by the measured M//i app ii ec j) 
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FIG. 2: (color online). The figure illustrates the evolution of 
the magnetization divided by the field strength, M(T,h)/h, 
with NLC order n and convergence towards experimental 
data. Panels a) and b) are for h = 0.2 T and h = 1.0 T, 
respectively, and with the field h applied along the [100] di- 
rection. Solid curves are for the NLC calculations and the 
number beside each curve corresponds to the order n up to 
which the calculations were carried to. Pluses are the demag- 
netization corrected experimental M(T, h) results divided by 
h. 



and N is the demagnetization factor, which depends on 
the sample geometry^ 



IV. NLC RESULTS: THERMODYNAMIC 
PROPERTIES 

The temperature dependence of the magnetization 
M(T, h) divided by the strength of the applied mag- 
netic field h along the cubic [100] direction and calcu- 
lated using the NLC method is shown in Fig. [5J Panels 
a) and b) of Fig. [2] show the results for a field of 0.2 
T and 1.0 T, respectively. The experimental values of 
M(T, h)/h, after demagnetization corrections, are shown 
for the same applied field values and are marked by black 
pluses. The number beside each curve labels the NLC 



order at which M(T,h) was calculated (see Section Hi) . 
NLC-0 (red curve, label "0") considers a single site clus- 
ter. It therefore does not incorporate the effect of the 
interactions in Eq.[T]between the Yb 3+ ions and thus cor- 
responds to the Yb 3+ single-ion property with g-tensor 
components g zz = 1.80 and g xy — 4.32. The significant 
differences between NLC-0 and the experimental data 
emphasizes the importance of interactions and concur- 
rent development of correlations below 10 K. The impor- 
tance of interactions on causing a departure of M(T, h)/h 
from its single-ion value below T < 10 K had previously 
been noted on the basis of a determination of the local 
susceptibility from polarized neutron measurement o 57 ' 58 
as well as a subsequent theoretical calculation of the local 
susceptibility^ using exchange parameters determined 
from RPA fits to the diffuse paramagnetic scattering^ 
However, none of these work o 57 ' 58 ' 66 had the quantitative 
descriptive power of the present NLC calculations. 

Compared to NLC-0, NLC-1 incorporates the effect of 
interactions and correlations at the scale of one tetra- 
hedron (see Fig. [I}. Focusing on the case h = 0.2 T, 
one observes little difference between NLC-0 and NLC-1 
above T ~ 2 K, but a large increase in M(T, h)/h in go- 
ing from NLC-1 to NLC-2 for T < 7 K. The results for 
NLC-1 requires the exact diagonalization of the Hamil- 
tonian over a single tetrahedron and therefore, incorpo- 
rates spatial spin-spin correlations extending only over a 
nearest- neighbor distance (see panel for n = 1 in Fig.QJ. 
In contrast, NLC-2 considers two tetrahedra (see panel 
for n = 2 in Fig. [IJ and therefore considers correlations, 
and associated fluctuations, reaching out to third near- 
est neighbors. As first noted by Ross et al.^ and further 
expanded upon by Applegate et al.^ the fluctuations of 
the joining spin/site of two tetrahedra mediates an effec- 
tive ferromagnetic third nearest-neighbor coupling ( J3 in 
Ref. [59[). This fluctuation-induced interaction promotes 
ferromagnetic correlations among the otherwise degener- 
ate classical 2-in/2-out spin ice states, an effect that we 
believe important to induce ferrimagnetic correlations in 
Yb2Ti 2 Q7i 20 ' 21 ' 59 The incorporation of this fluctuation 
process and induced effective J3 coupling only happens 
for NLC order n > 2 and is thus absent for NLC-1. 
We believe that the large increase of M(T, h)/h in going 
from NLC-1 to NLC-2 results from the ability of clus- 
ters of two-tetrahedra to support that type of fluctuation 
physics while a single tetrahedron cannot. One may want 
to argue that the observed large increase of the calculated 
M(T, h)/h, when going from NLC-1 to NLC-2, is further 
evidence that Yb2Ti 2 07, as described by the effective 
exchange parameters of Ref. [2(| , is either on its way or 
at the verge of developing spontaneous ferrimagnetic or- 
der at low-temperature. We return to this point below 
in Subsection [V] when we discuss the magnetic entropy, 
S(T, h), in nonzero magnetic field h. 

Considering the results of Fig. [2^), we observe that 
the NLC-2, 3 and 4 results are quite close to each other 
down to T ~ 1.8 K, the lowest temperature available 
for the experimental data. Most importantly, one should 
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note that NLC-2, 3 and 4, using the effective anisotropic 
exchange parameters {J e } of Ref. [13] already provide, 
without Euler summation, and without any adjustment 
of the parameters defining %qsi and Hz, a very good 
agreement with the experimental data. A similarly good 
agreement is also found between the experimental data 
and the NLC results for the highest NLC orders (n = 3, 4) 
for a field h = 1 T (see panel b) of Fig. |2J}. 

Having demonstrated that (i) interactions do play a 
strong effect in rcnormalizing M(T, h) for T < 10 K and 
that (ii) NLC orders n = 3 and n = 4 give a highly 
suitable description of the experimental data even for the 
weakest field /; = 0.2 T considered (see the discussion in 
Section HI)), we henceforth solely report results from the 
Euler transformation method (to order n = 3 and n = 4, 
see Eq. (|10p) to generate theoretical values of M(T, h)/h 
vs temperature and for various field h along the three 
high-symmetry cubic directions ([100], [110], and [111]). 

Figure [3] shows M(T,h)/h, obtained using the 4 th or- 
der Euler transformed NLC results, as a function of tem- 
perature T for fields h of strength 0.2, 1, 3, and 5 T 
oriented parallel to the [100], [110] and [111] directions. 
Again, the calculation employed the microscopic Hamil- 
tonian appropriate to Yb2Ti2C>7 as derived from inelastic 
neutron scattering data at high field 3^ Figure [3] shows 
4 th order NLC results for M (T, h) /h per Yb 3+ ion vs 
temperature T, with both quantities plotted on a log 
scale. The calculated magnetization is relatively inde- 
pendent of direction and levels off at low temperatures, 
with the temperature at which it levels off being lower 
with smaller applied magnetic fields. While the magne- 
tization is only weakly dependent on the direction of the 
applied magnetic field, such differences in the magneti- 
zation as a function of direction of field are only evident 
at low temperatures. 



V. MAGNETIZATION: CONFRONTING 
THEORY WITH EXPERIMENT 

Field dependent magnetization 

Figure U shows the experimentally-determined magne- 
tization as a function of temperature on a semi- log plot. 
Data is shown for applied magnetic field directions paral- 
lel to [100], [110], and [111], from top to bottom, and for 
field strengths of 0.2, 1, 3, and 5 T. Once again M(T, h) /h 
is normalized per Yb 3+ ion. Two sets of NLC expansion 
calculations using the Euler transformation are shown 
in Fig. 21 to order n = 3 are shown as the solid lines, 
and to order n — 4 are shown as the dot-dashed lines. 
Here again, the NLC expansions utilized the microscopic 
Hamiltonian derived from the high field inelastic neu- 
tron scattering (INS) data3£ with no adjustment in the 
parameters so determined. 

Several features are immediately clear from this com- 
parison between theory and experiment. First, the 
(Euler-transformed) NLC expansion results to both n — 
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FIG. 3: (color online). Calculated magnetizations as a func- 
tion of temperature for different field strengths and directions 
are shown using the 4 th order Euler transformed results from 
the NLC expansion (see. Eq. HSJ). The magnetic fields are 
applied along the [100], [110], and [111] directions. 



3 and n = 4 orders provide an excellent description of 
the magnetization for all field strengths and directions 
as well as all temperatures considered. The fact that 
a remarkable degree of quantitative agreement between 
theory and experiment is achieved without adjustment of 
the microscopic Hamiltonian appropriate to Yb2Ti207 
is strong validation of the determination of the micro- 
copic Hamiltonian^ As the INS measurements were per- 
formed at T = 0.03 K and applied magnetic fields of h = 
2 and 5 T with h parallel to [110], and the magnetization 
measurements are performed for T > 1.8 K and fields < 5 
T, we now have a very accurate description of Yb2Ti207 
using the same microscopic Hamiltonian over a remark- 
ably large region of its h — T phase diagram. 

Looking at the Euler-transformed NLC expansion re- 
sults for n = 3 (solid line) and n — 4 (dot-dashed line) in 
Fig. [4] one can see that, as expected, the two calculations 
are quite consistent with each other for T > 1.8 K, but 
depart from each other at lower temperatures. There is 
also a better agreement at lower temperatures between 
n = 3 and n = 4 NLC expansions for higher fields, in- 
dependent of the direction of the applied magnetic field, 
as anticipated on the basis of the general arguments pre- 
sented in Section |TTJ 



Other parametrizations of Hqsi 

There have been over the last four years, apart from 
Ref. [23], a number of other studies^ 3 -^^^^ combin- 
ing experiment and theory and which were aimed at de- 
termining the strength of the anisotropic interactions in 
Yb2Ti207. We believe that most, if not all, including 
Ref. [52| that was co-authored by one of us, are beset by 
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FIG. 4: (color online). Calculated magnetizations using the 
7i = 3 and n = 4 Euler-transformed M{T,h)/h NLC ex- 
pansion results (see Eq. (|10[)) are compared with the mea- 
sured magnetization versus temperature for fields applied 
along the [100], [ffO] and [111] directions. The parameters 
for H = "Hqsi + Hz are those determined form the spin wave 
results of Ref. [2^1 • Appropriate demagnetization corrections 
have been applied to all experimental data. Departures be- 
tween n — 3 and n = 4 NLC expansion results become noti- 
cible at temperatures < 1.8 K. 



significant drawbacks compared to the in-field inelastic 
neutron scattering measurements of Ross et al£Q In fact, 
one might have wondered whether the anisotropic ex- 
change parameters { J e } extracted in Ref. [2(| in magnetic 
fields of 5 Tesla might have suffered from large rcnormal- 
ization due to magnetoelastic effects. The good agree- 
ment between experimental and NLC results for zero 
magnetic field specific heat^ and the ones presented in 
this paper for the temperature and magnetic field de- 
pendent magnetization, M(t,h), provide compelling evi- 
dence that such renormalization, if it exists, is well within 
the experimental uncertainty of the estimated { J e } pa- 
rameters^ That being said, we now discuss each of the 
other works. 

Cao et al. used polarized neutron diffraction to extract 
the local susceptibility tensor of a number of RE2Ti2C>7 
pyrochlores, including Yb2Ti207, in an applied exter- 
nal field of 1 Tesla along the [100] crystallographic direc- 



tion 



57.58 



Their analysis was based on a mean-field theory 
that ignores the sublattice nature of the pyrochlore lat- 
tice and incorporates the effect of the local mean-field 
only via two independent coupling constants referred to 
as X z and Aj_. By construct, such an approach makes it 
difficult to reconstruct the microscopic exchange param- 
eters of the spin Hamiltonian. Furthermore, as they do 



not comment on this, it is not clear that their data anal- 
ysis took into account demagnetization effects. Malkin 
et al. reinvestigated the description of the bulk and local 
susceptibility of RE2Ti207 compounds, but now start- 
ing from a microscopic formulation, and incorporating 
an adequate sublattice structure in their model as well 
as including demagnetization corrections^ In the case 
of Yb2Ti2C>7, they reported being unable to describe the 
longitudinal site susceptibility \ 01 Cao e t oJ>- by us- 
ing a single set of crystal field parameters. It is per- 
haps important to note that Malkin et al. assumed that 
the bilinear anisotropic interactions between the mag- 
netic moment operators were symmetric, which amounts 
to neglecting Dzyaloshinskii-Moriya (DM) -like interac- 
tions. As was found in the work of Ross et al.^ the J z ± 
interactions, which originate from the DM interactions, 
are the second largest ones in Hqsi, and of almost the 
same magnitude as the largest J zz coupling. Hence, the 
neglect of the antisymmetric interactions with coupling 
J z ± for Yb2Ti2 07 is bound to cause difficulty in obtain- 
ing a quantitative description of the properties of this 
material. Consequently, because of the aforementioned 
caveats, the results obtained in Rcfs. 33ll57ll58j cannot 
not provide a quantitative account of the microscopic in- 
teractions at play in Yb2Ti 2 07. 

Thompson et al^ were the first to consider a micro- 
scopic theory which incorporates the single-ion crystal 
field Hamiltonian, H c f, for Yb 3+ in Yb2Ti207, the four 
symmetry-allowed bilinear nearest-neighbor interactions, 
Kf? between the Ji angular momentum operators along 
with the long-range magnetostatic dipole-dipole interac- 
tions^ In order to determine the K™, Thompson et al. 
used a random phase approximation (RPA) to calculate 
the diffuse (energy-integrated) neutron scattering inten- 
sity S(q) in the [hhl] scattering phase and compare with 
experimental measurements at a temperature T = 1.4 
With the Curie- Weiss temperature of Yb 2 Ti 2 7 
given by 9q w ~ 0.5 K ± 0.2 K, it would have seemed 
that RPA, which is in essence a mean-field scheme, may 
be a priori reasonably quantitatively accurate at a tem- 
perature of 1.4 K with a frustration/fluctuation level set 
by 0cw/T < 0.3. Forced by construction to describe 
short-range diffuse scattering, the RPA fit of Ref. [52[ 
provided a set of couplings K™ allowing for a good fit 
of S(q) and, by "retroactive consistency", a mean-field 
critical temperature T c ~ 1.1 K. When translating the 
determined Kf? parameters in the { J e } notation of Ross 
et al, one finds a large difference between the two sets. 
The authors of Ref. [39[ commented in the supplemen- 
tary material section of their paper that a possible rea- 
son for the failure of the model of Thompson et al^ to 



agree with the results of Ross et al^- is that the former 
work neglected multipolar interactions between the Ji 
operators beyond the considered bilinear ones. However, 
as discussed in the supplementary material of Ref. [52|, 
such critique 3 ^ would appear of little merit for the fol- 
lowing reason. Thanks to the essentially total isolation 
of the Kramers crystal field doublet of Yb 3+ from the 
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excited states, the ground crystal field doublet acts for 
an almost exact invertible unitary transformation of the 
Kij V Ji,uJj,v interactions, hence providing a one-to-one 
correspondence between the bilinear K"? couplings and 
the { Je] effective anisotropic exchange between effective 
spin-1/2 operators of Ross et al— _Ln other words, the 



bilinear AT"' 



?ives the J e SfS^ interac- 



Ji,uJj.v model of Ref. [52j can be viewed as 
a "high energy" model whose projection (with 'correct' 
values of the Kf- 1 couplings) 
tions in the model of Eq. (JTJ) [2C^ 

The difficulty with Thompson et al.'s results is readily 
understood on the basis of the mean-field theory results 
presented in Ref. [20|. The mean- field T™ f using the 
{ J e } parameters from Ref. [13] is approximately 3.5 K. 
Such a high T c mf compared to the 1.4 K RPA fit of S(q) 
means that the extracted A™ parameters suffer from a 
very significant and uncontrolled renormalization from 
thermal (and possibly quantum) fluctuations. Concerns 
that 1.4 K might be too low for a quantitative RPA fit of 
S(q) of Yb2Ti2C>7 and that fits at the higher temperature 
of 9.1 K data£2 might have been more appropriate had 
been expressed in Ref. [67|. Unfortunately, the diffuse 
signal at 9.1 K proved too weak to proceed. As demon- 
strated by Applcgatc et al. in their NLC calculations^ 
the anisotropic exchange parameters of Ref. [52| provide 
a poor description of the zero-field magnetic specific heat 
C(T) of Yb 2 Ti 2 7 . 

Chang et al. most recently reported their own estimate 
of the exchange parameters for Yb2Ti2 07<2 9 . They also 
employed an RPA scheme to fit the (polarized) neutron 
scattering intensity of the compound, but using an effec- 
tive spin-1/2 model rather than the bilinear A™ Ji. u Jj, v 
interactions plus H c { of Thompson et al£^ Perhaps be- 
lieving in the incorrectness of the latter description (see 
discussion two paragraphs above), they passed over in the 
body of their paper the opportunity to offer a critique of 
Thompson et al.'s model. That said, with the published 
evidence from Ref. [H that T c mf for Yb 2 Ti 2 7 may be 
as high as 3.5 K3i, and that the RPA fits of Thompson 
et al. at 1.4 K may thus be of questionable quantitative 
merit, it is surprising that Chang et al. did nevertheless 
proceed to use RPA to fit the polarized neutron scatter- 
ing at a temperature as low as 0.3 K, which is a mere 
25 percent higher than the experimental T c ~ 0.24 K. 
It is also not clear how the authors of Ref. [3!| find a 
similarity between their {J e } values and those of Ross et 
al., especially after evidence had been reported by Ap- 
plegate et aZ.^S that the { J e } of Ref. [39| fail to describe 
C (T) , while those of Ross et al. provide a quite adequate 
description of C(T)£^ Chang et al. propose a rational- 
ization of the low-temperature state of Yb2Ti207 on the 
basis of a Higgs-like phase and suggest positioning this 
material in a { J e } parameter space near a parent classical 
spin ice. Yet, at the same time, they overlook to com- 
ment on the inability of their microscopic model with 
its RPA-determined set of { J e } to describe C(T)£Z As 
in the case of Thompson et al.'s RPA analysis^ Chang 
et al.'s fit to the neutron scattering data provides, by 
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FIG. 5: (color online). The figure shows that M{T,h) from 
NLC calculations (orders n — 0, 1,4) using the {J e } parame- 
ters reported by Chang et a/.— fail dramatically to describe 
the experimental M(T,h) results (plus symbols). See text. 



the very consequence of using RPA at a temperature 
T <C T* nt , for a set of coupling parameters which are 
significantly and uncontrollably renormalized downward 
compared to the bare { J e }. This is illustrated in Fig. [5] 
where it is shown that NLC-1 and NLC-4 using the { J e } 
parameters of Chang et al. do not describe M(T,h) for 
temperatures T < 10 K, where correlations have barely 
started to develop^ In fact, down to T = 1 K, there is 
little difference between either NLC-1 and NLC-4 with 
the single-ion magnetization with all interactions turned 
off and given by the NLC-0 results. 

To conclude this discussion, it thus appears that only 
Ross et al.'s set of microscopic parameters^ consistently 
describe Yb2Ti207 for fields h < 5 Tcsla and tempera- 
tures T > 0.7 K. 



Field dependent specific heat 

While we are not aware of in-field specific heat, 
C{T,h), measurements on Yb2Ti 2 07 single crystals, we 
expect these to be soon carried out given the interest de- 
voted to this compound. One perspective as to why such 
measurements might be of interest is the following. The 
Er2Ti2C>7 pyrochlore antiferromagnet displays a transi- 
tion to long-range order at T c ~ 1.2 k£ i 69 ' 70 The appli- 
cation of a magnetic field along [110] ultimately destroys 
that order at T = 0, giving rise to a quantum phase 
transition at h c ~ 1.5 T. 69 ' 70 Specific heat measurements 
in non-zero field have been shown useful to characterize 
the evolution of this system in and out of the long-range 
ordered phase at h < 1.5 T and T < 1.2 K^ 9 .^ 

Measurements of C(T, h) in Yb 2 Ti 2 7 for T < 0.3 K 
may well be very interesting. This would be especially 
true if, once the sample dependence of T c has been un- 
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FIG. 6: (color online). The calculated molar specific heat us- 
ing NLC is shown as a function of temperature and magnetic 
field with field applied along, from top to bottom, the [100], 
[110] and [111] directions. These results were calculated us- 
ing the n = 3 and n — 4 Euler transformation NLC results 
using the {J e } parameters of Ref. [2(j. See the legends in top 
panel of Fig. [7] for values of the field strengths considered and 
legend for NLC order employed. 



derstood and T c 0.26 K indeed turns out to be a phase 
transition to a ferrimagnetically-ordered state, as sug- 
gested by some experiments and theory, but not all exper- 
iments i 35 ' 71 Thus, in anticipation of such measurements 
in Yb2Ti2C>7, we have used NLC expansion to calculate 
C(T,h). 

The specific heat, C(T,h)/R, calculated using the 
Euler-transformed NLC expansion to order n = 4 and the 
microscopic Hamiltonian determined for Yb2Ti207 from 
inelastic neutron scattering (INS)2fi is shown in Fig. [75] 
We present results for C(T,h)/R for a field parallel to 
[100], [110], and [111] from top to bottom, respectively. 
In all cases, calculations are presented for field strengths 
ranging from 0.2 T to 5 T. 

The trend is very similar for all three field directions. 
A broad peak is observed for a field strength h above ~ 
2 T, with the peak position moving to lower tempera- 
tures as h is decreased. While the NLC calculation is 
known to be inaccurate for temperatures less than ~ 1 
K^ 9 - and for appropriately low field strengths, the calcu- 
lated molar specific heat at the lowest applied magnetic 
fields is qualitatively consistent with that observed ex- 
perimentally in zero fieldj 3 -^ displaying a sharp anomaly 
at very low temperaures as well as a broad shoulder near 
2 K. 

Even though, at first sight, the dependence of ther- 



0.6 
0.4 
0.2 


-0.2 
0.6 
0.4 



— 0.2T 


NLC(E3) 


0.5T 


NLC(E4) 


LOT 




3.0T 




5.0T 






T(K) 

FIG. 7: (color online). The calculated molar entropy using 
NLC is shown as a function of temperature and magnetic field 
with field applied along, from top to bottom, the [100], [110] 
and [111] directions. These results were calculated using the 
n = 3 and n = 4 Euler transformation NLC results using the 
{J e } parameters of Ref. [2p( ] . 



modynamic properties looks independent of field direc- 
tion, there are subtle differences at low temperatures 
and fields, which could be taken as evidence of an ul- 
timate spontaneous ferrimagnetic ordering with M along 
[100] in zero field . 20 ' 21 i 39 ' 59 ' 71 On close inspection, one 
finds that for fields of 1 T and lower, C(T, h) is higher 
at temperatures around 1 K for the [100] direction and, 
as a result, the entropy S(T,h) is substantially lower. 
In fact, as shown in Fig. [7] the entropy along [100] ap- 
pears to become vanishingly small at temperatures of or- 
der 1 K, whereas there remains a residual entropy along 
the other two directions. The residual entropy remains 
largest along the [111] direction. This observation is con- 
sistent with the picture of a zero-field ordering in this 
system, characterized by M along one of the (100) crys- 
tallographic directions. 20 ' 21 ' 39 ' 59 In that case, applying a 
field along [100] selects one ordered state and the system 
merely exhibits a crossover to a paramagnetic state at a 
temperature of order 1 K. The entropy then displays an 
activated temperature dependence at low temperatures. 
In contrast, for small fields along [110] two degenerate 
ordered states remain and for small fields along [111], 
three degenerate states remain. Thus, there should con- 
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tinue to be a low temperature phase transition in suffi- 
ciently weak magnetic fields along either [110] or [111], 
with only below the transition the entropy going to zero. 
By its nature, NLC does not allow for our theoretical 
calculations to converge at low enough temperatures and 
magnetic fields to fully confirm this picture. It would 
thus be interesting to invesgitate its validity via C(T, h) 
experimental measurements. 

To summarize, the rather large anisotropic exchange 
interactions at play in Yb 2 Ti 2 07 lead to the develop- 
ment of a collective paramagnetic state at T < T™ f ~ 4 
K with 6qw ^™ f ; being thus somewhat "hidden" and 
difficult to quantitatively describe using standard text- 
book (e.g. RPA and other mean-field like) methods. In 
this material, and perhaps in other candidate quantum 
spin ice materials, the fit of bulk thermodynamic data 
such as C(T,h) and M(T,h) using NLC may provide a 
useful alternative. 



VI. DISCUSSION & CONCLUSION 

In this paper we have compared the results from nu- 
merical linked-cluster (NLC) calculations of the tem- 
parature and magnetic field dependent magnetization 
of Yb 2 Ti 2 07 with those obtained from experimental 
measurements on this material in a temperature range 
T e [1.8,20] K and magnetic field range h e [0.2,5] 
T. The NLC calculations were performed on a Hamil- 
tonian describing the interactions between pseudospins 
5=1/2 and characterized by effective anisotropic ex- 
change couplings determined via inelastic neutron scat- 
tering (INS) measurements in the polarized paramagnetic 
state of Yb 2 Ti 2 07i2£ The overall agreement between the 
NLC and the experimental results were found to be excel- 
lent (see Fig. SJ. Conversely, NLC magnetization results 
obtained using the exchange couplings determined from 
a random phase approximation analysis of the diffuse 
paramagnetic neutron scatterin g 39 - 52 ' 68 were found to be 
in significant disagreement with the experimental mag- 
netization measurements. The excellent argreement be- 
tween NLC and experimental measurements of M(T, h) 
shows that (i) the proposed nearest-neighbor exchange 
model 2 ^ is quantitatively accurate to describe Yb 2 Ti 2 07, 
that (ii) high-field inelastic scattering in the polarized 
paramagnetic state is a reliable method for extracting 
effective exchange parameters and that (iii) long-range 
dipolar interactions appear to not play an important role 
in the energetics of the system. The latter conclusion 
is perhaps a bit surprising given that long-range dipo- 
lar interactions play a key role in the physics of classi- 
cal (dipolar) spin ices . 30 ' 55 ' 56 One obvious difference is 
that, relative to the J zz effective exchange interaction, 
the nearest-neighbor contribution of the magnetostatic 
dipole-dipolc interactions in Yb 2 Ti 2 07 is approximately 
one order of magnitude smaller than in the classical dipo- 
lar spin ices . 11 ' 16 ' 30 ' 31 At the very least, one may con- 
sider that the nearest-neighbor contribution of the dipo- 



lar interactions are incorporated in the {J e } effective 
anisotropic exchange. Then, it thus appears that the 
perturbative long-range (beyond nearest-neighbor) part 
of the dipolar interactions plays no dramatic role in a 
field greater than 0.2 T and down to 1.8 K (this work) 
or, in zero field, down to approximately 0.7 K as found in 
Ref. [1^]. It may be that, if Yb 2 Ti 2 C>7 does possess a fer- 
rimagnetically ordered state with a magnetization along 
of the (100) axes, that the sole role of dipolar interactions 
is to weakly renormalize the critical temperature and the 
level of quantum fluctuations, along with inducing fer- 
rimagnetic domains. Yet, perhaps one should not be 
too expedient in assuming a generic irrelevance of long- 
range dipolar interactions of order 10 _1 compared to J zz 
for any candidate quantum spin ice material. One can 
imagine that dipolar interactions may play an important 
role, possibly inducing novel phases, in a material that 
would, based solely on its effective anisotropic exchange 
couplings {J e }, find itself at the boundary between the 
various semi-classical and intrinsically quantum phases 
identified in mean-field lattice gauge theories of quantum 
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spin icesi ' 

We believe that it would be, at this time, very inter- 
esting and most useful to carry out similar studies that 
combine inelastic neutron scattering and thermodynamic 
bulk measurements for other candidate quantum spin ice 
materials, with Yb 2 Sn 2 07 and Yb 2 Zr 2 07 being obvious 
choices. From the lessons learned in the present work, as 
well as from Ref. [5!|, we anticipate that NLC will con- 
tribute to developing a qualitative parametrization of mi- 
croscopic models of quantum spin ices and help pave the 
way for an ultimate understanding of these fascinating 
systems. However, for this program to be successful, the 
availability of large high resolution inelastic neutron scat- 
tering data sets will likely prove to be essential. In par- 
allel, in-field single crystal measurements which properly 
account for demagnetization effects will also be necessary. 
To end on Yb 2 Ti 2 07, we believe that, notwithstanding 
the evidence presented here and elsewher e) 37 ' 39 , it is by no 
means certai n 35 ' 71 at this time that long-range ferrimag- 
netic order develops below 0.26 K — in this very inter- 
esting material. More experiments on well-characterized 
samples are needed to settle this question. 
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